! subroutine to identify the bond list/type, angle list/type, dihedral
! list/type for a molecule based on its neighbor list and element name
subroutine identify()
use molecule
implicit none
   !----------------------------------------------------------------------------
   integer :: i, j, k, l, ni, nj, nk, nl, ip, jp, kp, lp, II, JJ, KK, LL
   integer :: idum1, idum2, ntype = 0, nmax
   !----------------------------------------------------------------------------
   !
   if ( allocated( types   ) ) deallocate( types    )
   if ( allocated( bonds   ) ) deallocate( bonds    )
   if ( allocated( angles  ) ) deallocate( angles   )
   if ( allocated( diheds  ) ) deallocate( diheds   )
   if ( allocated( typename) ) deallocate( typename )
   if ( allocated( bondname) ) deallocate( bondname )
   if ( allocated( anglename)) deallocate( anglename)
   if ( allocated( dihedname)) deallocate( dihedname)
   if ( allocated( charlist) ) deallocate( charlist )
   !
   ! Identify atom types by element name
   allocate( charlist(natom), types(natom) )
   ntype = 0
   do i = 1, natom
      types(i) = typeindex(element(i))
   enddo
   !
   natomtype = ntype
   allocate( typename(natomtype), typeid(natomtype) )
   forall (i=1:natomtype) typename(i) = trim(charlist(i))
   forall (i=1:natomtype) typeid(i) = i
   !
   ! Prepare for all other charlist
   nmax = natomtype * natomtype
   nmax = nmax * nmax
   deallocate( charlist )
   allocate( charlist( nmax ) )
   !
   ! bonds and bond types
   nbond  = count(neigh.gt.0)/2
   if ( nbond.gt.0 ) then
      allocate( bonds(3, nbond) )
      !
      ntype   = 0
      idum1   = 0
      do i = 1, natom
      do JJ = 1, count(neigh(:,i).gt.0)
         j  = neigh(JJ,i)
         if ( i.gt.j ) then
            if ( trim(element(i)).lt.trim(element(j)) ) then
               write(charflag, '(2A2)') element(i), element(j)
            else
               write(charflag, '(2A2)') element(j), element(i)
            endif
            idum1 = idum1 + 1
            bonds(:,idum1) = (/ i, j, typeindex(charflag) /)
         endif
      enddo
      enddo
      if ( idum1.ne.nbond) then
         write(*, '(I6, " bonds read, while ", I6, " bonds detected!")') nbond, idum1
         stop 'There must be something wrong with the neighbor list for bond!'
      endif 
      !
      nbondtype = ntype
      allocate( bondname(2,nbondtype), bondid(nbondtype) )
      do i = 1, nbondtype
         read(charlist(i), '(2A2)') bondname(1,i), bondname(2,i)
      enddo
      forall (i=1:nbondtype) bondid(i) = i
   endif
   !
   ! angles and angle types
   nangle = 0
   do i = 1, natom
      ni = count(neigh(:,i).gt.0)
      nangle = nangle + ni*(ni-1)/2
   enddo
   !
   if ( nangle.gt.0 ) then
      allocate( angles(4,nangle) )
      !
      ntype = 0
      idum1 = 0
      do i = 1, natom
         ni = count(neigh(:,i).gt.0)
         do JJ = 1, ni-1
         do KK = JJ+1, ni
            j  = neigh(JJ,i)
            k  = neigh(KK,i)
            if ( trim(element(j)).lt.trim(element(k)) ) then
               write(charflag, '(3A2)') element(j), element(i), element(k)
            else
               write(charflag, '(3A2)') element(k), element(i), element(j)
            endif
            idum1 = idum1 + 1
            angles(:,idum1) = (/ j, i, k, typeindex(charflag) /)
         enddo
         enddo
      enddo
      if ( idum1.ne.nangle) stop 'There must be something wrong with the neighbor list for angle!'
      !
      nangletype = ntype
      allocate( anglename(3,nangletype), angleid(nangletype) )
      do i = 1, nangletype
         read(charlist(i), '(3A2)') anglename(1,i), anglename(2,i), anglename(3,i)
      enddo
      forall (i=1:nangletype) angleid(i) = i
   endif
   !
   ! dihedrals and their types
   ndihed = 0
   do i = 1, natom
      ni = count(neigh(:,i).gt.0)
      do JJ = 1, ni
         j  = neigh(JJ,i)
         nj = count(neigh(:,j).gt.0)
         if ( j.gt.i ) cycle
         do KK = 1, ni
            k  = neigh(KK,i)
            if ( k.eq.j ) cycle
            do LL = 1, nj
               l  = neigh(LL,j)
               if (l.ne.i.and.l.ne.k) ndihed = ndihed + 1
            enddo
         enddo
      enddo
   enddo
   !
   if ( ndihed.gt.0 ) then
      allocate( diheds(5,ndihed) )
      !
      ntype = 0
      idum1 = 0
      !
      do i = 1, natom
         ni = count(neigh(:,i).gt.0)
         do JJ = 1, ni
            j  = neigh(JJ,i)
            nj = count(neigh(:,j).gt.0)
            if ( j.gt.i ) cycle
            do KK = 1, ni
               k  = neigh(KK,i)
               if ( k.eq.j ) cycle
               do LL = 1, nj
                  l  = neigh(LL,j)
                  if (l.ne.i.and.l.ne.k) then
                     if ( trim(element(i)).eq.trim(element(j)) ) then
                        if ( trim(element(k)).lt.trim(element(l)) ) then
                           write(charflag, '(4A2)') element(k), element(i), element(j), element(l)
                        else
                           write(charflag, '(4A2)') element(l), element(j), element(i), element(k)
                        endif
                     else if ( trim(element(i)).lt.trim(element(j)) ) then
                        write(charflag, '(4A2)') element(k), element(i), element(j), element(l)
                     else
                        write(charflag, '(4A2)') element(l), element(j), element(i), element(k)
                     endif
                     idum1 = idum1 + 1
                     diheds(:,idum1) = (/ k, i, j, l, typeindex(charflag) /)
                  endif
               enddo
            enddo
         enddo
      enddo
      ndihedtype = ntype
      allocate( dihedname(4,ndihedtype), dihedid(ndihedtype) )
      do i = 1, ndihedtype
         read(charlist(i), '(4A2)') dihedname(1:4,i)
      enddo
      forall (i=1:ndihedtype) dihedid(i) = i
   endif
            
contains
   !
   integer function typeindex(flag)
   use molecule, only: charlist
   implicit none
      !-------------------------------------------------------------------------
      character (len=*), intent(in) :: flag
      integer                       :: it
      !-------------------------------------------------------------------------
      typeindex = 0
      do it = 1, ntype
         if ( trim(flag).eq.trim(charlist(it)) ) then
            typeindex = it
            return
         endif
      enddo
      !
      if ( typeindex.eq.0 ) then
         ntype = ntype + 1
         typeindex = ntype
         charlist(ntype) = trim(flag)
      endif
   return
   end function typeindex
end subroutine
